All runs Ulva & Hypnea rETRmax Analysis, Script Chunks, and Plots
This is the analysis of all the Ulva lactuca salinity and nutrient experiments conducted on the lanai in St. John 616 from September 2021 to October 2022. These experiments incorporated four paired salinity and nutrient treatments with three temperatures. Each of the first four runs produced an n = 2 and was repeated initially 8 times for a total of n = 16. Data gaps were identified and filled in February, April, and October 2022. This output reflects all data totaling five treatments for Ulva lactuca and six for Hypnea musciformis. Update 2/18/23: all rETRmax are from the Y.II>0.1 modification in dataset
Packages loaded:
library(lme4)
library(lmerTest)
library(afex)
library(effects)
library(car)
library(MuMIn)
library (dplyr)
library(emmeans)
library(DHARMa)
library(performance)
library(patchwork)
library(rstatix)
#for plots and tables
library(ggplot2)
library(ggpubr)
library(forcats)
library(RColorBrewer)
library(tidyverse)
library(sjPlot)
library(sjmisc)
library(mmtable2)
library(gt)
library(purrr)
library(stringr)
library(tidyr)
Load this dataset normalized to quantum efficiency of photosynthesis per Silsbe and Kromkamp 2012
all_runs_photosyn_data <- read.csv("../data_input/hyp_ulva_all_runs_ek_alpha_normalized.csv")
Assign several variables as factors
all_runs_photosyn_data$Run <- as.factor(all_runs_photosyn_data$Run)
all_runs_photosyn_data$Temperature <- as.factor(all_runs_photosyn_data$Temp...C.)
all_runs_photosyn_data$Treatment <- as.factor(as.character(all_runs_photosyn_data$Treatment))
all_runs_photosyn_data$deltaNPQ <- as.factor(all_runs_photosyn_data$deltaNPQ)
Subset the species for output. Use Day 9 for final analysis. Remove treatment 2.5 from Ulva dataset
ulva <- subset(all_runs_photosyn_data, Species == "ul" & RLC.Day == 9 & Treatment != 2.5)
ulva$treatment_graph[ulva$Treatment == 0] <- "1) 35ppt/0.5umol"
ulva$treatment_graph[ulva$Treatment == 1] <- "2) 35ppt/14umol"
ulva$treatment_graph[ulva$Treatment == 2] <- "3) 28ppt/27umol"
ulva$treatment_graph[ulva$Treatment == 3] <- "5) 18ppt/53umol"
ulva$treatment_graph[ulva$Treatment == 4] <- "6) 11ppt/80umol"
ulva$treatment_graph[ulva$Treatment == 2.5] <- "4) 28ppt/53umol"
#ULVA rETRmax________________________________________________________________
Make a histogram of the data
ulva %>% ggplot(aes(rETRmaxYpoint1)) +
geom_histogram(binwidth=5, fill = "#5BB300", color = "black", size = 0.25, alpha = 0.85) +
theme_bw()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
Run model without interaction between the treatments and temperature
ulva_retrmax_model <- lmer(formula = rETRmaxYpoint1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva)
Make residual plots of the data for ulva
hist(resid(ulva_retrmax_model))
plot(resid(ulva_retrmax_model) ~ fitted(ulva_retrmax_model))
qqnorm(resid(ulva_retrmax_model))
qqline(resid(ulva_retrmax_model))
Check the performance of the model, make a table, plot effects, and get r2
performance ::check_model(ulva_retrmax_model)
## Variable `Component` is not in your data frame :/
r.squaredGLMM(ulva_retrmax_model)
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
## R2m R2c
## [1,] 0.5908239 0.6867304
summary(ulva_retrmax_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## rETRmaxYpoint1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +
## (1 | RLC.Order)
## Data: ulva
##
## REML criterion at convergence: 1828.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8318 -0.6181 0.0792 0.5191 3.3964
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plant.ID (Intercept) 22.232 4.715
## Run (Intercept) 6.655 2.580
## RLC.Order (Intercept) 5.189 2.278
## Residual 111.306 10.550
## Number of obs: 240, groups: Plant.ID, 95; Run, 7; RLC.Order, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 34.850 2.978 9.421 11.702 6.36e-07 ***
## Treatment1 20.925 3.209 7.309 6.521 0.000272 ***
## Treatment2 21.808 3.209 7.309 6.796 0.000209 ***
## Treatment3 38.133 3.209 7.309 11.883 4.83e-06 ***
## Treatment4 39.456 3.209 7.309 12.296 3.80e-06 ***
## Temperature27 -4.434 2.511 30.820 -1.766 0.087337 .
## Temperature30 -2.246 2.296 66.339 -0.978 0.331519
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trtmn1 Trtmn2 Trtmn3 Trtmn4 Tmpr27
## Treatment1 -0.640
## Treatment2 -0.640 0.775
## Treatment3 -0.640 0.775 0.775
## Treatment4 -0.640 0.775 0.775 0.775
## Temperatr27 -0.403 0.003 0.003 0.003 0.003
## Temperatr30 -0.391 0.000 0.000 0.000 0.000 0.475
tab_model(ulva_retrmax_model, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)
| rETRmaxYpoint1 | ||||||
|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | Statistic | p | df |
| (Intercept) | 34.85 | 2.98 | 28.98 – 40.72 | 11.70 | <0.001 | 229.00 |
| Treatment [1] | 20.93 | 3.21 | 14.60 – 27.25 | 6.52 | <0.001 | 229.00 |
| Treatment [2] | 21.81 | 3.21 | 15.49 – 28.13 | 6.80 | <0.001 | 229.00 |
| Treatment [3] | 38.13 | 3.21 | 31.81 – 44.46 | 11.88 | <0.001 | 229.00 |
| Treatment [4] | 39.46 | 3.21 | 33.13 – 45.78 | 12.30 | <0.001 | 229.00 |
| Temperature [27] | -4.43 | 2.51 | -9.38 – 0.51 | -1.77 | 0.079 | 229.00 |
| Temperature [30] | -2.25 | 2.30 | -6.77 – 2.28 | -0.98 | 0.329 | 229.00 |
| Random Effects | ||||||
| σ2 | 111.31 | |||||
| τ00 Plant.ID | 22.23 | |||||
| τ00 Run | 6.66 | |||||
| τ00 RLC.Order | 5.19 | |||||
| ICC | 0.23 | |||||
| N Run | 7 | |||||
| N Plant.ID | 95 | |||||
| N RLC.Order | 6 | |||||
| Observations | 240 | |||||
| Marginal R2 / Conditional R2 | 0.591 / 0.687 | |||||
plot(allEffects(ulva_retrmax_model))
Construct null model to perform likelihood ratio test REML must be FALSE
ulva_retrmax_treatment_null <- lmer(formula = rETRmaxYpoint1 ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
ulva_retrmax_model2 <- lmer(formula = rETRmaxYpoint1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
anova(ulva_retrmax_treatment_null, ulva_retrmax_model2)
## Data: ulva
## Models:
## ulva_retrmax_treatment_null: rETRmaxYpoint1 ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## ulva_retrmax_model2: rETRmaxYpoint1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## npar AIC BIC logLik deviance Chisq Df
## ulva_retrmax_treatment_null 7 1986.1 2010.4 -986.04 1972.1
## ulva_retrmax_model2 11 1873.0 1911.3 -925.52 1851.0 121.04 4
## Pr(>Chisq)
## ulva_retrmax_treatment_null
## ulva_retrmax_model2 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ulva_retrmax_temperature_null <- lmer(formula = rETRmaxYpoint1 ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
ulva_retrmax_model3 <- lmer(formula = rETRmaxYpoint1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
anova(ulva_retrmax_temperature_null, ulva_retrmax_model3)
## Data: ulva
## Models:
## ulva_retrmax_temperature_null: rETRmaxYpoint1 ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## ulva_retrmax_model3: rETRmaxYpoint1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## npar AIC BIC logLik deviance Chisq Df
## ulva_retrmax_temperature_null 9 1871.9 1903.2 -926.95 1853.9
## ulva_retrmax_model3 11 1873.0 1911.3 -925.52 1851.0 2.8515 2
## Pr(>Chisq)
## ulva_retrmax_temperature_null
## ulva_retrmax_model3 0.2403
Plots
ulva %>% ggplot(aes(treatment_graph, rETRmaxYpoint1)) +
geom_boxplot(size=0.5) +
geom_point(alpha = 0.75, size = 3, aes(color = Temperature), show.legend = TRUE) +
labs(x="salinity/nitrate", y= "rETRmax (μmols electrons m-2 s-1)", title= "A", subtitle = "Ulva lactuca") +
scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "18ppt/53umolN", "11ppt/80umolN")) +
ylim(-1, 170) + stat_mean() +
scale_color_manual(values = c("#295102", "#7CB950", "#BDE269")) +
geom_hline(yintercept=0, color = "red", size = 0.5, alpha = 0.5) +
theme_bw() +
theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05),
plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))
Summarize the means for rETRmax
ulva %>% group_by(Treatment) %>% summarise_at(vars(rETRmaxYpoint1), list(mean = mean))
## # A tibble: 5 × 2
## Treatment mean
## <fct> <dbl>
## 1 0 32.6
## 2 1 53.2
## 3 2 54.1
## 4 3 70.4
## 5 4 71.8
ulva %>% group_by(Run) %>% summarise_at(vars(rETRmaxYpoint1), list(mean = mean))
## # A tibble: 7 × 2
## Run mean
## <fct> <dbl>
## 1 1 64.7
## 2 2 58.4
## 3 3 64.2
## 4 3.5 66.3
## 5 4 61.2
## 6 6 35.8
## 7 6.5 29.4
ulva %>% group_by(Treatment, RLC.Day) %>% summarise_at(vars(rETRmaxYpoint1), list(mean = mean))
## # A tibble: 5 × 3
## # Groups: Treatment [5]
## Treatment RLC.Day mean
## <fct> <int> <dbl>
## 1 0 9 32.6
## 2 1 9 53.2
## 3 2 9 54.1
## 4 3 9 70.4
## 5 4 9 71.8
Add growth rate to the dataset and subset by species
growth_rate <- read.csv("/Users/Angela/src/work/limu/algal_growth_photosynthesis/data_input/all_runs_growth_011723.csv")
growth_rate$Species <- as.factor(growth_rate$Species)
growth_rate$treatment <- as.factor(growth_rate$treatment)
growth_rate$lunar.phase <- as.factor(growth_rate$lunar.phase)
Make a new column for weight change (difference final from initial)
growth_rate$growth_rate_percent <- (growth_rate$final.weight - growth_rate$Initial.weight) / growth_rate$Initial.weight * 100
gr_ulva <- subset(growth_rate, Species == "Ul" & treatment != 2.5)
ulva$growth_rate <- round((gr_ulva$final.weight - gr_ulva$Initial.weight) / gr_ulva$Initial.weight * 100, digits = 2)
ulva$lunar.phase <- (gr_ulva$lunar.phase)
Plot a regression between the photosynthetic independent variables of interest and growth rate
ulva_growth_etr_graph <- ggplot(ulva, aes(x=rETRmaxYpoint1, y=growth_rate)) +
geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = Treatment)) +
geom_smooth(method = "lm", col = "black") + theme_bw() +
labs(title = "Ulva lactuca rETRmax vs Growth Rate", x = "rETRmax (μmols electrons m-2 s-1)",
y = "growth rate (%)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor()
ulva_growth_etr_graph
## `geom_smooth()` using formula = 'y ~ x'
HYPNEA____________________________________________________________________________________________________________
There was no D9 RLC for hm6-4 on 11/12/21 but had to remove hm6-4 from 10/9/21 below to match growth data because replicate appeared dead
hypnea <- subset(all_runs_photosyn_data, Species == "hm" & RLC.Day == 9 & uid != "2021-10-09_hm6-4")
hypnea$treatment_graph[hypnea$Treatment == 0] <- "1) 35ppt/0.5umol"
hypnea$treatment_graph[hypnea$Treatment == 1] <- "2) 35ppt/14umol"
hypnea$treatment_graph[hypnea$Treatment == 2] <- "3) 28ppt/27umol"
hypnea$treatment_graph[hypnea$Treatment == 3] <- "5) 18ppt/53umol"
hypnea$treatment_graph[hypnea$Treatment == 4] <- "6) 11ppt/80umol"
hypnea$treatment_graph[hypnea$Treatment == 2.5] <- "4) 28ppt/53umol"
For Hypnea remove hm6-4 on 11/12 that had no d9 RLC (final weight 0.1017) and hm6-4 on 10/9/21 because it was white and also white and dead hm1-1 on 10/12/22 and hm1-2 on 4/29/22 causing issues of influential observations hm1-2 for both rETRmax and Ek – leaving them in dataset because no good reason to believe not good data
gr_hypnea <- subset(growth_rate, Species == "Hm" & final.weight != 0.1017 & growth_rate_percent > -87.96837)
hypnea$growth_rate <- round((gr_hypnea$final.weight - gr_hypnea$Initial.weight) / gr_hypnea$Initial.weight * 100, digits = 2)
Make a histogram for hypnea
hist(hypnea$rETRmaxYpoint1, main = paste("Hypnea musciformis"), col = "maroon", labels = TRUE)
hypnea %>% ggplot(aes(rETRmaxYpoint1)) +
geom_histogram(binwidth=5, fill = "#7D0033", color = "black", size = 0.25, alpha = 0.85) +
theme_bw()
Run model for rETRmax
hyp_retrmax_model <- lmer(formula = rETRmaxYpoint1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea)
Make residual plots of the data
hist(resid(hyp_retrmax_model))
plot(resid(hyp_retrmax_model) ~ fitted(hyp_retrmax_model))
qqnorm(resid(hyp_retrmax_model))
qqline(resid(hyp_retrmax_model))
Check the performance of the model, get r2, make table, and plot effects
performance::check_model(hyp_retrmax_model)
## Variable `Component` is not in your data frame :/
r.squaredGLMM(hyp_retrmax_model)
## R2m R2c
## [1,] 0.2320523 0.6199556
summary(hyp_retrmax_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## rETRmaxYpoint1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +
## (1 | RLC.Order)
## Data: hypnea
##
## REML criterion at convergence: 2223
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.16586 -0.53708 -0.07062 0.46478 3.15105
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plant.ID (Intercept) 40.71 6.381
## Run (Intercept) 63.05 7.941
## RLC.Order (Intercept) 13.08 3.617
## Residual 114.48 10.700
## Number of obs: 286, groups: Plant.ID, 143; Run, 8; RLC.Order, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 56.646 6.278 6.289 9.023 7.93e-05 ***
## Treatment1 1.574 7.125 5.373 0.221 0.833
## Treatment2 -2.910 7.125 5.373 -0.408 0.699
## Treatment2.5 10.625 10.052 4.639 1.057 0.342
## Treatment3 2.476 7.125 5.373 0.348 0.741
## Treatment4 4.733 7.133 5.398 0.663 0.534
## Temperature27 -14.339 2.938 30.237 -4.880 3.21e-05 ***
## Temperature30 -15.962 2.569 46.255 -6.214 1.35e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trtmn1 Trtmn2 Trt2.5 Trtmn3 Trtmn4 Tmpr27
## Treatment1 -0.777
## Treatment2 -0.777 0.953
## Treatmnt2.5 -0.551 0.485 0.485
## Treatment3 -0.777 0.953 0.953 0.485
## Treatment4 -0.776 0.952 0.952 0.485 0.952
## Temperatr27 -0.220 0.001 0.001 0.000 0.001 0.000
## Temperatr30 -0.209 0.000 0.000 0.000 0.000 0.003 0.468
tab_model(hyp_retrmax_model, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)
| rETRmaxYpoint1 | ||||||
|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | Statistic | p | df |
| (Intercept) | 56.65 | 6.28 | 44.29 – 69.01 | 9.02 | <0.001 | 274.00 |
| Treatment [1] | 1.57 | 7.13 | -12.45 – 15.60 | 0.22 | 0.825 | 274.00 |
| Treatment [2] | -2.91 | 7.13 | -16.94 – 11.12 | -0.41 | 0.683 | 274.00 |
| Treatment [2.5] | 10.62 | 10.05 | -9.16 – 30.41 | 1.06 | 0.291 | 274.00 |
| Treatment [3] | 2.48 | 7.13 | -11.55 – 16.50 | 0.35 | 0.728 | 274.00 |
| Treatment [4] | 4.73 | 7.13 | -9.31 – 18.78 | 0.66 | 0.508 | 274.00 |
| Temperature [27] | -14.34 | 2.94 | -20.12 – -8.55 | -4.88 | <0.001 | 274.00 |
| Temperature [30] | -15.96 | 2.57 | -21.02 – -10.91 | -6.21 | <0.001 | 274.00 |
| Random Effects | ||||||
| σ2 | 114.48 | |||||
| τ00 Plant.ID | 40.71 | |||||
| τ00 Run | 63.05 | |||||
| τ00 RLC.Order | 13.08 | |||||
| ICC | 0.51 | |||||
| N Run | 8 | |||||
| N Plant.ID | 143 | |||||
| N RLC.Order | 6 | |||||
| Observations | 286 | |||||
| Marginal R2 / Conditional R2 | 0.232 / 0.620 | |||||
plot(allEffects(hyp_retrmax_model))
Construct null model to perform likelihood ratio test REML must be FALSE
hypnea_retrmax_treatment_null <- lmer(formula = rETRmaxYpoint1 ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
hypnea_retrmax_model2 <- lmer(formula = rETRmaxYpoint1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
anova(hypnea_retrmax_treatment_null, hypnea_retrmax_model2)
## Data: hypnea
## Models:
## hypnea_retrmax_treatment_null: rETRmaxYpoint1 ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## hypnea_retrmax_model2: rETRmaxYpoint1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## npar AIC BIC logLik deviance Chisq Df
## hypnea_retrmax_treatment_null 7 2283.4 2309.0 -1134.7 2269.4
## hypnea_retrmax_model2 12 2279.3 2323.2 -1127.7 2255.3 14.039 5
## Pr(>Chisq)
## hypnea_retrmax_treatment_null
## hypnea_retrmax_model2 0.01536 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hypnea_retrmax_temperature_null <- lmer(formula = rETRmaxYpoint1 ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
hypnea_retrmax_model3 <- lmer(formula = rETRmaxYpoint1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
anova(hypnea_retrmax_temperature_null, hypnea_retrmax_model3)
## Data: hypnea
## Models:
## hypnea_retrmax_temperature_null: rETRmaxYpoint1 ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## hypnea_retrmax_model3: rETRmaxYpoint1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## npar AIC BIC logLik deviance Chisq Df
## hypnea_retrmax_temperature_null 10 2307.6 2344.1 -1143.8 2287.6
## hypnea_retrmax_model3 12 2279.3 2323.2 -1127.7 2255.3 32.248 2
## Pr(>Chisq)
## hypnea_retrmax_temperature_null
## hypnea_retrmax_model3 9.941e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Plots
hypnea %>% ggplot(aes(treatment_graph, rETRmaxYpoint1)) +
geom_boxplot(size=0.5) +
geom_point(alpha = 0.75, size = 3, aes(color = Temperature), show.legend = TRUE) +
labs(x="Treatment", y= "rETRmax (μmols electrons m-2 s-1)", title= "B", subtitle = "Hypnea musciformis") +
scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) +
ylim(-1, 170) + stat_mean() +
scale_color_manual(values = c("#9C0627", "#BB589F", "#F4B4E2")) +
geom_hline(yintercept=0, color = "red", size = 0.5, alpha = 0.5) +
theme_bw() +
theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05),
plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))
Plot a regression between the photosynthetic independent variables of interest and growth rate
hypnea_growth_etr_graph <- ggplot(hypnea, aes(x=rETRmaxYpoint1, y=growth_rate)) +
geom_point(alpha = 0.5, size = 3, aes(color = Treatment), show.legend = TRUE) +
geom_smooth(method = "lm", col = "black") + theme_bw() +
labs(title = "Hypnea musciformis rETRmax vs Growth Rate", x = "rETRmax (μmols electrons m-2 s-1)",
y = "growth rate (%)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor(label.y = 150)
hypnea_growth_etr_graph
## `geom_smooth()` using formula = 'y ~ x'
Summarize the means
hypnea %>% group_by(Treatment) %>% summarise_at(vars(rETRmaxYpoint1), list(mean = mean))
## # A tibble: 6 × 2
## Treatment mean
## <fct> <dbl>
## 1 0 46.5
## 2 1 49.2
## 3 2 44.7
## 4 2.5 57.2
## 5 3 50.1
## 6 4 52.6
hypnea %>% group_by(Treatment, RLC.Day) %>% summarise_at(vars(rETRmaxYpoint1), list(mean = mean))
## # A tibble: 6 × 3
## # Groups: Treatment [6]
## Treatment RLC.Day mean
## <fct> <int> <dbl>
## 1 0 9 46.5
## 2 1 9 49.2
## 3 2 9 44.7
## 4 2.5 9 57.2
## 5 3 9 50.1
## 6 4 9 52.6